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INTRODUCTION  AND  OBJECTIVES 


Practical  applications  of  computational  fluid  dynamics  within  the  foreseeable 
future  must  necessarily  utilize  some  form  of  turbulence  modeling.  Computations 
based  on  the  Reynolds-averaged  equations  of  motion  require  all  turbulence  trans¬ 
port  of  momentum  and  energy  to  be  modeled;  while  large-eddy  simulations  require 
the  subgrid-scale  turbulence  to  be  modeled.  In  both  eases,  the  uncertainties  in 
modeling  turbulence  within  the  viscous  sublayer  constitute  a  major  weak  link  in 
the  overall  numerical  computation. 

During  the  past  two  decades  a  great  deal  of  new  experimental  information  has 
been  assembled  on  the  physics  of  organized  eddy  structures  in  turbulent  flow,  espe¬ 
cially  within  the  viscous  sublayer  (see,  for  example,  the  recent  review  of  Cantwell, 
1081).  Yet  it  has  not  been  possible  thus  far  to  incorporate  this  body  of  physical 
information  within  the  framework  of  Reynolds-averaged  turbulence  modeling.  The 
reason  is  fundamental,  reflecting  a  well-known  limitation  in  the  Reynolds-average 
approach  which  begins  by  time  averaging  the  dynamic  equations  of  motion.  In  this 
initial  mathematical  step  important  physical  aspects  of  organized  eddy  structures, 
such  as  phase  relationships  and  coherent  structure  dynamics,  are  obliterated  irre¬ 
versibly.  Consequently,  some  totally  different  approach  is  required  if  the  observed 
physics  of  coherent  eddy  motions  are  to  be  incorporated  within  the  framework  of  a 
turbulence  model. 

Quite  recently  a  new  approach  has  been  developed  for  modeling  viscous  sub¬ 
layer  turbulence  in  incompressible  flow  without  heat  transfer  (Chapman  and  Kuhn 
(1981).  Their  method  models  directly  the  essential  organized  eddy  structures  ob- 

i 

served  in  experiments.  The  principal  steps  in  this  ‘coherent-structure  modelling* 
are:  first,  to  model  velocity  boundary  conditions  at  the  outer  edge  of  the  viscous  sub¬ 
layer,  then  to  compute  time-dependent  dynamics,  and  finally  to  time  average  com- 


puted  results.  Thus,  time  averaging  is  the  last  operation  performed  on  computed 
dynamics,  rather  than  the  first  operation  performed  on  the  dynamic  equations.  This 
initial  effort,  although  far  from  fully  developed,  already  has  been  surprisingly  suc¬ 
cessful  in  modeling  some  important  characteristics  of  viscous  sublayer  turbulence. 

The  over  all  objective  of  the  present  research  is  to  develop,  using  the  Navier- 
Stokes  equations,  a  computational  model  of  viscous  sublayer  turbulence  applicable 
to  flow  with  heat  transfer  and  to  compressible  flow.  Specific  tasks  within  this 
objective  are  to  utilize  the  model  to  compute  the  distribution  of  turbulent  Prandtl 
number  within  the  viscous  sublayer  (a)  for  fluids  of  various  molecular  Prandtl 
number,  and  (b)  for  flows  with  adverse,  zero,  and  favorable  streamwise  pressure 
gradients. 

Navier-Stokes  computations  of  the  turbulent  Prandtl  number  Prt  are  significant 
because  experimental  techniques  have  been  unable  to  provide  reliable  measurements 
within  the  viscous  sublayer.  This  parameter,  of  course,  affects  turbulent  heat 
transfer.  Since  by  definition, 

Prt- 

an  experimental  determination  would  require  accurate  measurements  near  a  wall 
of  the  Reynolds  stress  on,  the  correlation  Tv  between  temperature  #  and  normal 
velocity,  the  mean  velocity  gradient  dTJ/dy,  and  the  mean  temperature  gradient 
dT/dy.  To  date,  it  has  not  been  possible  to  make  reliable  measurements  of  these  four 
quantities  in  the  viscous  sublayer.  The  band  of  uncertainty  from  one  set  of  data  for 
air  (Simpson  et  al.  (1970)),  is  illustrated  in  fig.  1.  Other  sets  of  experimental  data, 
e.g.  Fulachier  (1972),  indicate  that  the  overall  uncertainty  band  is  even  broader  than 
this  one  set  of  measurements  would  indicate.  Experiments  are  unable  to  determine 
whether  the  values  of  Prt  near  a  wall  are  large,  small,  or  intermediate.  Experiments 
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also  have  not  been  able  to  define  how  Prt  varies  with  molecular  Pr  or  with  pressure 
gradient.  We  believe  that  Navier-Stokes  computations  can  shed  much  light  on  these 
uncertainties. 

MOTIVATIONS  FOR  COHERENT  STRUCTURE  MODELING  OF 
VISCOUS  SUBLAYER  TURBULENCE 

Coherent-structure  modeling  departs  markedly  from  Reynolds-averaging  model¬ 
ing,  and  thereby  offers  the  potential  of  some  entirely  new  advances  in  turbulence 
computation.  Both  the  motivations  and  the  payoff  for  this  type  of  research  are 
quite  different  from  those  for  more  conventional  turbulence  modeling.  In  addition 
to  the  specific  motivations  outlined  above  for  the  objectives  of  the  present  research, 
there  are  other  important  motivations  that  will  affect  different  areas  of  turbulence 
computation.  It  is  appropriate,  therefore,  to  outline  some  of  these. 

Three  other  motivations  for  developing  a  realistic  coherent-structure  model  of 
viscous  sublayer  turbulence  are: 

1.  -  To  Strengthen  a  major  weak  link  in  present  Reynolds-Average  closure  Schemes. 

Conventional  turbulence  closure  methods  (e.g.  k—t,  other  2-equation  methods, 
and  Reynolds-Stress  transport  methods)  all  employ  differential  transport  equa¬ 
tions  that  are  “modeled”  forms  of  exact  transport  equations.  The  exact  equa¬ 
tion  for  the  transport  of  dissipation,  for  example,  contains  very  complex  tur¬ 
bulence  terms  involving  pressure-velocity  correlations,  triple  correlations  of 
velocity  gradients,  correlations  of  second  derivatives  of  velocity,  and  correla¬ 
tions  between  pressure  gradient  and  velocity  gradient.  Such  correlations  are  not 
measurable  with  present  experimental  technology.  Because  modeled  equations 
for  free  turbulence  yield  demonstrably  incorrect  results  near  a  wall,  various  ad 
hoc  functions  (up  to  5  in  number)  are  added  in  an  effort  to  mend  this  shortcom- 


ing.  Without  any  guide  or  test  basis  from  experiment,  the  inevitable  conse¬ 
quence  is  that  different  models  with  different  ad  hoc  functions  yield  different 
results  for  technically  important  quantities  such  as  skin  friction  in  Sows  with 
pressure  gradient  (Patel  et.  al  1981).  If,  however,  a  realiatic  Navier-Stokes 
computational  model  were  developed  for  the  time  dependent  dynamics  of  vis¬ 
cous  sublayer  flow,  then  all  the  terms  in  the  exact  transport  equations  could 
be  computed.  When  compared  with  corresponding  modeled  terms,  a  rational 
formulation  would  be  obtained  of  the  wall  damping  functions  for  integrating 
turbulence  models  across  the  viscous  sublayer.  This  would  strengthen  a  major 
weak  link  in  conventional  turbulence  models.  Reynolds-average  modeling,  of 
course,  will  be  the  mainstay  of  practical  turbulent  flow  computation  for  years 
to  come. 

2.  -  To  provide  a  simple  test  flow  against  which  subgrid  scale  models  of  turbulence 
in  large  eddy  simulations  can  be  tested.  The  simple  flow  of  homogeneous  shear, 
for  example,  computed  in  detail  from  the  time-dependent  Navier-Stokes  equa¬ 
tions  (Rogallo  1977)  has  been  used  effectively  as  a  test  flow  for  assessing  sub¬ 
grid  scale  turbulence  models  for  application  in  the  outer  region  of  a  turbulent 
boundary  layer  or  in  free  shear  layers  (e.g.,  Clark  et.  al.  1979,  McMillan 
and  Ferziger  1979,  Shirani  1981).  A  realistic  coherent-structure  model  of  vis¬ 
cous  sublayer  turbulence,  likewise  computed  from  the  time-dependent  Navier- 
Stokes  equations,  could  similarly  be  used  to  test  subgrid-scale  models  for  large 
eddy  simulations  of  turbulent  flow  in  the  region  adjacent  to  a  wall. 

3.  -  To  provide  a  guide  for  modeling  the  lower  boundary  conditions  for  the  outer 
turbulent  region  in  large  eddy  simulations  of  boundary  layer  flow.  If  a  realistic 
boundary  condition  of  this  type  could  be  developed,  the  entire  viscous  sublayer 
could  be  modeled  rather  than  directly  computed.  The  required  computer  power 
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for  large  eddy  simulations  at  high  Reynolds  numbers  would  thereby  be  reduced 
by  a  very  large  factor  (about  10*;  Chapman,  1980). 

In  summary  of  the  motivations  for  developing  coherent-structure  turbulence 
models,  it  is  clear  that  the  potential  payoffs  are  significant  and  of  much  broader 
scope  than  the  specific  objectives  of  the  present  research.  Benefits  are  anticipated 
to  conventional  Reynolds-average  models  of  near-term  application,  as  well  as  to 
future  large  eddy  simulation  models  of  longer  range  application. 


FORMULATION  OF  COMPUTATIONAL  MODEL 
Differential  Equations 

The  numerical  scheme  selected  to  solve  the  Navier-Stokes  equations  is  an 
adaptation  of  the  Pulliam  and  Steger  (1980)  code  for  compressible  flow.  Velocities 
in  the  physical  z,  y,  ~i  domain  are  u,  v,  w.  In  the  computational  (,  ij,  (,  domain 
the  dependent  variables  are  p,  pu,  pv,  pW,  and  y[e,  +  (u2  +  02  +  w2)/ 2],  where  e,-  is 
the  internal  energy  per  unit  mass.  Pulliam  and  Steger  use  z  as  the  coordinate 
normal  to  the  surface,  and  y  as  the  spanwise  coordinate.  In  order  to  minimize  code 
changes,  we  use  this  coordinate  system  in  the  numerical  computations  and  then,  at 
the  end,  change  J  to  y,w  to  v,y  to  z,  and  v  to  w  in  order  to  present  results  in  the 
more  familiar  boundary  layer  coordinate  designations.  We  employ  simple  cartesian 
coordinates  in  the  streamwise  z-  and  spanwise  p-directfcros,  and  a  stretched  mesh 
in  the  normal  7  direction  in  order  to  concentrate  points  both  near  the  wall  and  the 
outer  edge.  The  metrics  for  the  transformation  are  —  1,  qy  ■»  1/Af+,  j*  —  1/A7+ 
where  Ap+  and  A7+  are  the  mesh  spacings  in  wall  units  in  the  spanwise  and  normal 
directions  respectively. 

The  same  basic  physical  approximation  is  made  as  in  Chapman  and  Kuhn 
(1981),  namely,  that  coherent  sublayer  eddy  structures  are  highly  elongate  stream- 
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wise.  In  our  case  we  assume-in  accordance  with  experimental  observations  (e.g. 
Iritani  et.  al  (1981))  —  that  the  temperature  eddies  are  also  highly  elongate  stream- 
wise.  In  the  notation  of  Pulliam  and  Steger,  the  corresponding  divergence  form  of 
the  Navier-Stokes  equations  is 


dtq  -  d„{Fv  —  F)  +  d<{Gv  -  G)  +  d^F,  (1) 

where  the  vectors  F,  Fv,  G,  Gv  and  q  are  identical  to  the  conventional  ones  used 
by  Pulliam  and  Steger  for  the  case  of  orthogonal  coordinates  corresponding  to  the 
metrics  of  equation  (1).  However,  the  vector 


Fp  -  J"1 


0 

Pe 

0 

0 

-^iPe&We. 


(2) 


differs  from  theirs  due  to  our  assumption  of  highly  elongate  eddies  in  the  streamwise 
direction.  Thus,  streamwise  derivatives  other  than  the  mean  streamwise  pressure 
gradient  dpe/dx  are  neglected  relative  to  derivatives  in  the  spanwise  and  normal 
directions. 

This  basic  approximation,  which  involves  computation  of  3  velocity  components 
in  2  space  directions,  has  been  variously  termed  “2.5D  flow”,  or  “slender  turbulence 
theory”,  or  “slender  eddy  theory.” 

Since  Pulliam  and  Steger  used  the  thin-layer  form  of  the  Navier-Stokes  equa¬ 
tions,  we  write,  in  order  to  retain  a  tridiagonal  structure, 


d(Fv-F) 

dq 

d(Gv  —  G) 


dP 

dq 

ad 

dq 


Ac  +  ib 

Be  +  Bg 


(3) 
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where  the  matrices  A#  and  B0  contain  no  cross  derivative  terms,  while  At  and  Bx 
contain  all  such  terms.  The  detailed  structure  of  these  matrices  is  given  in  Appendix 
I.  Our  basic  factored  algorithm  becomes 


(/-  |W(BS)A4"+I  =  lh[6t,Fn  +  dfG"]  +  \hdnA**r 

+  +  -  A  qn 


(4) 


this  differs  from  the  Pulliam-Steger  algorithm  by  the  addition  of  the  last  four 
terms  on  the  right  hand  side.  Of  these,  the  last  one  corresponds  to  a  three-point 
backward  Euler  difference  scheme  for  time  derivatives  (added  for  accuracy  and 
stability  reasons  explained  later);  the  next  to  last  term  corresponds  to  body  pressure 
gradients  imposed  by  the  large  scale  eddies,  while  the  third  and  fourth  from  last 
represent  the  additional  terms  required  for  the  full  Navier-Stokes  equations  in  place 
of  the  thin-layer  equations  used  by  Pulliam  and  Steger. 

BOUNDARY  CONDITIONS 

The  primary  element  in  formulating  our  computational  model  involves  the 
construction  of  realistic  time-  and  space-dependent  boundary  conditions  for  the 
fluctuating  temperature  f  and  the  fluctuating  velocity  components  u,  v,  v  at  the 
outer  edge  of  the  viscous  sublayer.  Various  experimental  observations  are  used  as  a 
guide  to  this  construction.  Computations  using  the  time-dependent  Navier-Stokes 
equations  are  then  made  over  a  length  of  time  long  enough  for  the  various  statistical 
quantities,  e.g.,  mean  values,  correlation  coefficients,  rms  fluctuating  intensities, 
etc.,  to  reach  a  periodic  state  within  the  viscous  sublayer. 

Our  general  approach  is  an  extension  of  that  employed  by  Chapman  and  Kuhn 
(1081)  for  incompressible  flow  without  heat  transfer.  They  considered  two  coherent 
components  of  turbulent  motion  at  the  outer  edge  of  the  viscous  sublayer  (lower 
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edge  of  the  logarithmic  region).  One  component  represents  the  small  scale  eddies 
(SSE)  responsible  for  the  principal  production  of  turbulence  and  Reynolds  stress.  A 
second  component  represents  the  organized  large  scale  eddies  (LSE)  which  interact 
with  the  small  scale  eddies.  They  used  a  simple  initial  construction  for  purposes 
of  illustrating  the  potential  and  main  characteristics  of  this  type  of  turbulence 
modeling.  With  all  quantities  expressed  in  conventional  dimensionless  wall  variables, 
their  velocity  boundary  conditions  for  the  outer  edge  of  the  viscous  sublayer  were 

Component  1-SSE  Component  2-LSE 


ue+  —  2oj  sin  N\T  sinf  +  ^2(o2  —  of  )sin(NM2T  +  ^*2) 

ve+  =  —2/1  sin  NiT  sin  f 

w«+  20cosN1T cos f  +  \J T1  —  fP)nvlr~T  +  ^*,2) 


where  a,  /?,  -j  are  the  rms  intensities  of  fluctuation  of  ue+,  ve,  «()  respectively,  N\  is 
the  mean  frequency  of  SSE  burst  events  (ejection/sweep  events),  Nu 2  is  the  mean 
frequency  of  LSE,  f  ss  2nx/\  is  a  dimensionless  spanwise  distance  variable,  T  is  the 
dimensionless  time,  X  is  the  mean  spacing  between  high-speed  or  low-speed  streaks, 
«i/o  is  the  uu  correlation  coefficient  (0.45),  and  d«2  and  2  are  phase  angles.  The 
flow  model  is  periodic  both  in  time  and  span.  Hence  the  boundary  conditions  for 
the  spanwise  side  walls  of  the  computational  domain  were 

««+(*,  0,0  — Ue+(»,X,«) 

0, 0  —  »r+(»,  X,  t) 

««+(»,  0,  *)  —  ««+(»,  X,<) 
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and  at  the  surface 


Ue+(0,  z,  t)  =  ve+(0,  z,  t)  =  we+(0,  z,  t)  =  0 

This  construction  of  velocity  boundary  conditions  was  guided  by  experimental 
observations  of  organized  sublayer  structure.  For  example,  the  sin  factor  for  ue+ 
corresponds  to  the  observation  of  high-  and  low-speed  streaks  spaced  spanwise  a 
mean  distance  of  X  apart,  while  the  factors  sinj,  sinj,  cosf  in  the  ue+,  t>e+,  *>e+ 
equations,  respectively,  correspond  to  a  simple  modeling  of  the  observations  of 
contra-rotating  vortical  motion  near  a  wall.  The  absence  of  a  phase  angle  in  time 
between  ue+  and  ve+  corresponds  to  the  observation  that  the  u  —  v  correlation  in 
the  log  region  is  maximum  with  zero  time  delay  between  their  respective  signals. 
The  18CP  phase  difference  between  ue+  and  ve+  corresponds  to  observations  from 
conditional  samples  in  the  sublayer  that  ti  and  v  are  18(f  out  of  phase  during  the 
Reynolds-stress  intensive  ejection-sweep  event.  The  9Cf  phase  difference  in  time 
between  ve  and  we  corresponds  to  the  observation  that  the  derivative  (dv2/dy+)e 
at  the  outer  edge  of  the  viscous  sublayer  is  zero.  A  comparison  of  a  number  of 
turbulence  characteristics  computed  from  this  simple  model  with  those  measured 
in  experiments  showed  surprisingly  good  agreement.  This  model  was  developed  for 
flows  with  a  mean  streamwise  pressure  gradient,  although  numerical  computations 
were  made  by  Chapman  and  Kuhn  only  for  the  relatively  small  pressure  gradients 
that  exist  in  incompressible  pipe  and  channel  flow. 

The  above  two-component  velocity  boundary  conditions  at  the  outer  edge 
of  the  viscous  sublayer  are  fully  compatible  with  the  concept  of  “active”  and 
“inactive”  components  of  turbulent  motion  in  the  log  region  as  characterized  by 
Townsend  (1961)  and  Bradshaw  (1967).  Component  1  is  active,  involving  small 
scale  eddies,  producing  the  Reynolds  stress,  being  rotational,  and  being  dependent 
on  wall  variables.  Component  2  is  inactive,  involving  large-scale  eddies,  producing 


energy  but  no  Reynolds  stress,  being  irrotational,  and  being  dependent  on  outer 
variables.  There  is  much  experimental  information  that  turbulent  flow  in  the  log 
region,  and  hence  at  the  outer  edge  of  the  viscous  sublayer,  comprises  these  two 
distinct  types  of  motion. 

The  computational  model  of  Chapman  and  Kuhn  is  not  regarded  as  sufficiently 
realistic  for  purposes  such  as  outlined  in  the  section  above  (Hi  motivations.  Their 
results  do  illustrate,  however,  that  the  general  approach  employed  has  the  potential 
of  eventually  leading  to  a  sufficiently  realistic  model.  Towards  this  end,  it  is  believed 
that  at  least  three  additional  aspects  of  turbulence  physics  will  have  to  be  introduced 
into  the  model.  First,  production  of  Reynolds  stress  should  be  intermittent  in 
time,  rather  than  sinusoidal- like  as  in  boundary  conditions  (5).  Second,  conditional 
samples  for  v(t)  in  ejection/sweep  events  indicate  that  the  v2  turbulent  energy 
contained  in  these  SSE  events  should  be  only  a  fraction  of  the  total  v 2  energy, 
rather  than  the  entire  amount  as  in  equation  (5).  Third,  visualizations  of  viscous 
sublayer  flow  exhibit  a  mixture  of  order  and  disorder,  rather  than  coherent  order 
only.  To  date,  in  work  supported  by  another  contract,  progress  has  been  made  on 
modifying  the  outer  edge  velocity  boundary  conditions  to  account  for  the  first  two 
of  these  three  physical  aspects.  We  will  incorporate  this  progress  into  our  velocity 
boundary  conditions.  In  addition,  for  the  heat  transfer  case  we  must  also  construct 
an  appropriate  boundary  condition  for  the  fluctuating  temperature  9e  at  the  outer 
edge  of  the  viscous  sublayer. 

One  important  guide  to  the  construction  of  boundary  conditions  is  provided 
by  measurements  of  spectral  density  /.  Fortunately,  spectra  for  f,  u,  v,  and  v,  have 
been  obtained  by  Fulachier  (1972)  at  y-positions  near  the  outer  edge  of  the  viscous 
sublayer.  Data  for  ye+  <=»  40  (interpolated  between  his  measurements  at  y+  =  31 
and  77)  are  shown  in  figure  2.  Since 


-  10  - 


kfd(lu  k)=  l 


roo  r oo 

l 

the  areas  under  curves  of  u2kfu,  v2kfV)  and  w2kfw  versus  logit  are  proportional  to 
the  relative  amounts  of  kinetic  energy  in  these  velocity  components.  Along  the  wave 
number  scale  k,  pips  are  shown  corresponding  to  LSE  eddy  scales  of  j t/6  (where  6 
is  the  boundary  layer  thickness)  and  to  SSE  eddy  scales  of  jt/X.  It  is  seen  that, 
whereas  both  u  and  8  exhibit  eddy  scales  ranging  from  large  (0(6))  to  small  (0(X)], 
v  exhibit  scales  ranging  only  from  medium  [0  (several  X)]  to  small.  Our  boundary 
conditions  will  be  constructed  to  be  as  consistent  as  possible  with  experimental  data 
of  this  type. 

A  relatively  simple  method  (developed  under  a  separate  contract)  is  used  to 
simulate  the  intermittent  character  of  Reynolds-stress  production.  Since  Sv  is 
produced  only  by  the  active  component  of  small  scale  eddies  (SSE),  only  the  time 
functions  sin  and  cmN\T  in  equation  (1)  need  be  modified.  For  this  purpose 
we  use  a  simple  Fourier  series  approximation  H(j>)  to  a  rectangular  pulse  function, 
where  <t>  =  N\T. 


For  M  terms  in  the  approximating  Fourier  series,  we  have 

M 

H(t)  =  ^  (an  cos  +  bn  sin  n^) 

1 

B5=|E(*2.+*5) 

an  ™  — -  sin  nXi  —  sin  sXj 

X2 

6„  =  —  ^(cos  nXi  —  1)  +  (cos  nX\  —  1)1 

jrnLA2  J 

values  of  M  in  the  range  of  3  to  5  are  anticipated.  Since  the  skewness  of  the 
rectangular  pulse  function  is 

c  =  (X2-Xj) 

*  VXiXtJXi+Xi 

we  can  build  skewness  into  simply  by  having  Xi  different  from  X\.  Initial 
computations,  however,  are  planned  for  zero  skewness  at  the  outer  edge,  and  with 
X\  =  X2  *=  0.25«-,  since  hot-wire  measurements  have  indicated  that  the  principal 
uv  is  produced  in  about  25%  of  the  time,  with  relative  St?  quiescence  the  remaining 
75%  of  time. 

In  conventional  wall  variables,  the  outer  edge  boundary  conditions  constructed 
to  date  include  terms  representing  small  scale  eddies  (SSE)  of  length  scale  X,  large 
scale  eddies  (LSE)  of  length  scale  >  10X,  and  medium  scale  eddies  (MSE)  of  inter¬ 
mediate  length  scale  3X.  Computations  are  made  within  a  domain  covering  3  span- 
wise  cells  each  X  in  width.  The  LSE  are  treated  as  being  of  sufficiently  large  scale 
that  their  time-dependent  velocity  component  does  not  vary  across  the  spanwise 
extent  of  the  computational  domain.  The  equation  of  continuity  precludes  having 
a  LSE  component  in  the  boundary  condition  for  ve. 
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SSE 


LSE 

Scale  >10X 


MSE 
Scale  3X 


Scalfi-X 


«e  =  >/2ai,F0(^)sinf  +  ^2(a2  -  a^)s\n(Nu2T  +  ^o2) 

vt  =  y/2^\Fv(4>)s\a  f  +  2^2  sin  (3N.2r)  sin  | 

—  +  y4(72  -  £?)sin  ^  j 

(8) 

tfe  =  \/2ai,Fj(^)sinf  +  ^2(o2  -  of)sin(N.2T  +  tn) 


where,  in  addition  to  the  quantities  already  defined,  a  ss 
in  the  large-eddy  component  of  temperature,  and 


is  a  phase  angle 


F,(t)  =  -Fu(*)  =  Fv(<t>)  = 
*w\9)  , - -  = 

]/hH+  +  +»  1) 


H(4>) 

\f&&) 


are  the  intermittent  time  functions  normalized  such  that  F5  =  1.  The  boundary 
conditions  for  the  spanwise  side  walls  of  the  comutational  domain  are 


«e+(v,0,  t)  =  «e+(y,3X,f) 
«*+(y»o.t)*=w*+(if.3X,0 
»«+(*,  0,  t)  *=■  v«+(y,  3X,  I) 


Thus  the  computation  allows  for  flow  interaction  between  the  center  cell  and  the 
two  outside  cells  unconstrained  by  specific  boundary  conditions. 
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As  in  the  model  of  Chapman  and  Kuhn:  a,  fi,  a  are  determined  from 
experimental  data  on  the  intensities  of  turbulence  fluctuations  in  wall  units  at  the 
outer  edge  of  the  viscous  sublayer  (a  «  2,^  l,i  »  !J,s  «  1.3);  4«,2»  4«2> 

are  determined  by  computer  trial  to  conform  to  the  mean-velocity  law  of  the  wall, 
mean-temperature  law  of  the  wall,  and  u-velocity  skewness,  respectively;  4mi  is 
determined  analytically  by  the  requirement  that  (dv2/dy)e  is  zero;  aj  is  determined 
by  the  correlation  coefficient  0.45  =  ~(RUv)e  =  aifa/aft;  oy  is  determined  by  the 
correlation  coefficient  0.8  -(Rus)e  =  +  ^^(ft2  - 

and  0x1  f)  is  determined  from  conditional  samples  of  the  ejection/sweep  event.  This 
latter  determination  is  made  by  assuming  that  the  peak  to  peak  amplitude  ratio 

<  v  >PPa  /  <  u  >PPa,  from  the  conditional  samples  is  eqnal  to  the  ratio  fiifai, 
and  then  using  the  above  equation  for  (Auti)e  to  compute  both  py/fi  and  ay/a .  An 
independent  determination  of  ay/a  has  been  made  from  direct  measurements  of 
the  fraction  of  total  u2  energy  that  exists  during  bursts,  and  hence  also  of  P\fp 
from  the  equation  for  (RUv)c  Still  a  different  determination  has  been  made  from 

<  ®  >ppa  /  <  «  >ppa—  Pi/<*i  coupled  with  the  equation  for  Results  of  these 

determinations  (made  on  a  separate  contract)  are: 


Data  Source 

Method 

Pl/P 

Chen  &  Blackwelder  (1978) 

<  v  >  1  <  u  >  plus  Ruv  «=  —.45 

.53 

Nakagawa  &  Nezu  (1981) 

<  v  >  /  <  u  >  plus  RvV  «—  —.45 

.72 

Blackwelder  &  Kaplan  (1979) 

<  v  >  /  <  u  >  plus  RnV  —  —.45 

.49 

Kim  (1982) 

<  u  >  /  <  v  >  from  LES  computations, 
plus  Ruv  =  -.45 

Kim,  Kline  &  Reynolds 

Fraction  of  total  u2  energy  during 

.59 

(1971) 

bursts  (.68),  plus  RnV  -.45 

Blackwelder  &  Kaplan  (1979) 

<«>/<«>  plus  —  -.45 

.64 

On  average,  the  value  of  Py/P  appears  to  be  near  0.6. 

It  is  to  be  noted  that  the  structure  of  boundary  conditions  (8)  comprises  in  equal 
proportions  both  ejection/sweep  events  (where  a  sweep  follows  an  ejection),  and 
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sweep/ejection  events  (where  an  ejection  follows  a  sweep).  The  recent  experiments 
of  Johansson  and  Alfredsson  (1982)  clearly  reveal  that  both  types  of  event  occur. 
Their  experimental  VITA  data  for  the  smallest  threshold  values  and  the  longest 
integration  times  employed,  indicate  these  two  event  types  to  be  about  equally 
numerous.  Both  types  have  also  been  observed  visually  by  Offen  and  Kline  (1974) 
who  use  the  term  "cleansing  sweep”  in  an  ejection/sweep  event  to  distinguish  from 
sweeps  in  a  sweep/ejection  event. 


STATUS  of  research 

As  outlined  above,  research  on  formation  of  the  differential  equations  and  the 
initial  trial  boundary  conditions  is  completed.  We  also  have  (1)  essentially  completed 
a  survey  of  various  existing  theories  for  turbulent  Prandtl  number,  and  have  (2) 
made  most  of  the  programming  modifications  necessary  to  utilize  the  Pulliam- 
Steger  code  in  our  computational  model.  Some  details  of  progress  in  these  latter 
two  areas  are  presented  in  the  sections  which  follow. 

Survey  of  various  theories  for  Prt 

Without  a  reliable  guide  from  experiment  as  to  how  Prt  should  be  modeled  near 
a  wall,  the  various  theoretical  models  advanced  to  date  have  differed  greatly.  This  is 
illustrated  in  figures  3,  4,  and  5  showing  our  compilation  of  various  theories  for  fluids 
of  molecular  Prandtl  number  0.72  (air),  6  (water),  and  1000  (oil),  respectively.  Our 
survey  of  existing  theories  is  near  completion,  although  not  yet  fully  complete.  The 
various  theories  explored  clearly  demonstrate  that  a  very  wide  range  of  uncertainty 
exists  for  Prt  in  the  viscous  sublayer.  Present  theories  do  not  attempt  to  define  how 
Prt  may  vary  with  pressure  gradient. 
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Programming  of  Navier-Stokes  Code 

Our  plan  is  to  develop  a  single  Navier-Stokes  code  which  can  be  used  either  for 
low-speed  subsonic  flows  (M  <  <  1)  with  heat  transfer,  or  for  high  Mach  number 
flows.  Since  the  Pulliam-Steger  code  was  written  for  the  thin-layer  approximation 
to  the  Navier-Stokes  equations,  a  number  of  programming  modifications  are  re¬ 
quired  to  adapt  this  code  to  the  computation  of  viscous  sublayer  dynamics.  These 
re: 


1.  Change  variables  from  free  stream  reference  to  nondimensional  wall  variables. 

2.  Add  subroutine  for  the  boundary  conditions  on  velocity  and  temperature  fluctuations 
at  the  outer  edge  of  the  viscous  sublayer. 

3.  Change  the  time-derivative  algorithm  to  a  stable,  second-order  accurate  algo¬ 
rithm. 

4.  Add  the  spanwise  viscous  stress  terms,  plus  the  associated  cross  derivatives 
terms. 

5.  Add  a  global  streamwise  pressure  gradient  term,  associated  with  the  large  scale 
eddies,  to  the  right  hand  side  of  the  basic  algorithm. 

6.  Add  a  block  periodic  solver  for  the  spanwise  direction. 

7.  Add  appropriate  nonperiodic  terms  in  the  spanwise  boundary  conditions  for  p, 

P,  T. 

An  outline  of  the  progress  status  of  each  of  these  seven  items  is  given  in  the  sections 

which  follow. 

1.  Completed 

2.  Subroutine  coding  of  the  boundary  conditions  for  the  initial  test  case  has  been 
completed.  This  case  comprises  the  same  velocity  boundary  conditions  as 
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used  by  Chapman  and  Kuhn  (1081),  plus  a  temperature  boundary  condition 
representing  small-scale  and  large-scale  eddy  components.  The  subroutine  will 
be  used  to  check  results  against  the  previous  results  of  Chapman  and  Kuhn 
for  low-speed  subsonic  flow  without  heat  transfer.  Once  a  check  is  obtained,  a 
modified  boundary-conditions  subroutine  will  be  written  to  represent  boundary 
conditions  (8)  which  simulate  intermittent  Reynolds  stress  production,  rather 
than  the  sinusoidal-like  production  contained  in  the  Chapman-Kuhn  boundary 
conditions  (5). 

3.  Completed.  In  order  to  test  the  algorithm  accuracy  for  the  time  derivative 
terms,  an  exact  analytical  solution  was  derived  for  incompressible  oscillating 
shear  flow  with  heat  transfer.  This  solution  is  an  extension  of  the  Chapman- 
Kuhn  solution  for  flow  without  heat  transfer.  Mathematical  details  of  the 
analytical  solution  are  presented  in  Appendix  II.  It  was  found  that  neither  of 
the  two  option  algorithms  in  the  Pulliam-Steger  code  (i.e.,  Euler  implicit  first 
order,  or  trapezoidal  second  order)  were  sufficiently  accurate  and  stable.  With 
helpful  guidance  from  Dr.  Pulliam,  a  three-point  backward  Euler  algorithm 
was  coded  and  found  to  be  satisfactory.  The  very  close  agreement  obtained 
with  this  second-order  accurate  algorithm,  between  numerical  computations 
and  the  exact  analytical  solution,  is  illustrated  in  figures  6  and  7  for  the 
adiabatic-wall  case  and  the  heat-transfer  case,  respectively.  Since  oscillating 
shear  flow  with  heat  transfer  involves  viscous  and  heat-conduction  terms  as 
well  as  time-derivative  terms,  this  test  flow  provides  a  check  on  the  accuracy 
of  the  numerics  for  all  three  types  of  terms.  It  does  not  provide,  however,  a 
check  on  the  accuracy  of  the  convection  terms  which  are  zero  in  oscillating 
shear  flow. 

4.  Completed 

5.  Completed 
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6.  Completed 

7.  Not  yet  undertaken.  This  coding  modification  apparently  is  not  needed  for  low 
speed  subsonic  flows.  The  necessary  code  changes  are  not  planned  to  be  un¬ 
dertaken  until  we  begin  to  model  high  Mach  number  flows.  With  fluctuations 
in  spanwise  body  pressure  gradient,  all  three  velocity  components  remain  peri¬ 
odic  spanwise,  but  the  state  variables  p,  p,  0  do  not.  The  linear  body  pressure 
gradients  create  step  discontinuities  at  the  spanwise  boundary  points.  This  can 
induce  numerical  oscillations  near  the  boundaries  with  the  central  differencing 
scheme  used.  At  low  subsonic  Mach  numbers  these  discontinuities  are  very 
small,  so  that  the  block  periodic  solver  can  be  used  with  little  problem.  At 
high  Mach  numbers,  however,  the  relative  pressure  fluctuations  Ap/p  increase, 
and  the  magnitude  of  the  step  discontinuities  can  increase  sufficiently  to  cause 
numerical  problems  if  periodic  boundary  conditions  are  imposed  on  the  state 
variables.  Our  plan  to  handle  this  for  high  speed  flows  is  to  decompose  p,  p  and 
T  into  a  spanwise  periodic  component  determined  by  small  scale  eddies,  plus  a 
nonperiodic  component  determined  by  large  scale  eddies.  The  latter  would  be 
included  as  a  right-hand-side  term,  and  the  overall  set  of  equations  cast  into 
a  form  for  which  periodic  boundary  conditions  are  applicable. 
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PROGRESS  SUMMARY 


During  the  first  contract  year  (CY  1982)  we  have  made  progress  as  follow* 
Survey  of  theories  for  Pr»  essentially  complete. 

Exact  analytical  solution  derived  for  incompressible  oscillating  shear  Sow  with 
heat  transfer- 

Using  analytical  solution,  satisfactory  test  conducted  on  accuracy  of  numerical 
algorithms  for  time-dependent,  heat-conduction,  and  viscous  terms. 

Most  necessary  programming  modifications  to  the  Pulliam-Steger  Navier-Stokes 
code  completed  (but  not  yet  debugged). 

A  boundary  condition  for  fluctuating  temperature  has  been  formulated  for  initial 
trial  computations. 

Some  of  the  subsequent  steps  to  be  conducted  in  CY  1988  are: 

Complete  necessary  programming  modifications  to  Pulliam-Steger  code  (e.g.  various 
statistical  space-time  averages  for  fluctuation  intensities,  correlations,  skewness, 
flatness,  dissipation,  Reynolds  stresses,  etc. 

Debug  code  and  run  check  case  for  incompressible  flow  without  heat  transfer 
using  same  boundary  conditions  as  employed  by  Chapman  and  Kuhn. 

Make  computer  runs  to  develop  temperature  boundary  conditions  that  are  as 
realistic  as  possible,  testing  computations  against  experimental  data. 

Make  computer  runs  for  fluids  of  various  molecular  Prandtl  number  and  stream- 
wise  pressure  gradient. 
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APPENDIX  I 


DECOMPOSITION  MATRICES 

The  derivatives  dF/dg  and  dGjdg  are  each  decomposed  into  two  matrices,  one 
(subscript  zero)  that  does  not  contain  cross  derivative  terms,  and  another  (subscript 
x )  that  contains  all  such  terms. 


-r - A,  +  Ax 

-  Bo  +  BX 
dq 


where  the  matrix  Aa  is 
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and  the  matrix  Bx  is 
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APPENDIX  H 


An  exact  analytical  solution  has  been  derived  for  oscillating  shear  flow  with  heat 
transfer  in  an  incompressible  fluid  with  constant  properties.  Since  the  momentum 
equations  are  independent  of  the  thermal  energy  equation  in  incompressible  flow,  the 
velocity  field  is  the  same  as  for  the  case  without  heat  transfer  derived  by  Chapman 
and  Kuhn  (1981).  With  y  as  the  coordinate  normal  to  the  infinite  flat  surface, 
and  u  and  v  as  the  velocity  components  in  the  streamwise  and  spanwise  directions, 
respectively,  we  have 

u  =  A^s'in  nt  —  t~kn!l  sin(nt  —  itny)J  +  cy 
w  =  £?[sin(mt  +  (^)  —  e~kmV  sin(mf  —  kmy) | 


where 


The  oscillating  frequencies  are  n  in  the  streamwise  direction  and  m  in  the  spanwise 
direction,  and  the  mean  velocity  profile  is  B  =  cy.  For  this  velocity  field,  the 
solution  to  the  thermal  energy  equation 
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Figure  1.  -  Envelope  of  experimental  uncertainty  for  turbulent  Prandtl  number 
from  data  of  Simpson  et.  al.  (1970). 


Figure  2.  Spectral  density  of  the  three  components  of  velocity  fluctuation  at 
y+  ■  40  from  data  of  Fulachier  (1972). 
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Figure  3.  Some  theoretical  models  of  turbulent  Prandtl  number 
for  a  molecular  Prandtl  number  of  0.72  (air). 
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Figure  4.  Some  theoretical  models  of  turbulent  Prandtl  number 
for  a  molecular  Prandtl  number  of  6  (water) . 
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Figure  5.  Some  theoreClcal  models  of  turbulent  Prandtl  number 
for  a  molecular  Prandtl  number  of  1000  (oil) . 
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Figure  6.  Comparison  of  numerical  computations 
with  exact  analytical  results  for  os¬ 
cillating  shear  flow;  Adiabatic  wall. 


